The Idea of this is to emulate the information gathered from nibabel and spm_dicom_orient.py into one handy jupyter notebook with the explanations along with the sympy output.
This page describes the DICOM orientaion. The sections will be as follows:
First we define the standard DICOM patient-based coordinate system. This is what DICOM means by $x$, $y$ and $z$ axes in its orientation specification. From section C.7.6.2.1.1 of the From section C.7.6.2.1.1 of the DICOM object definitions (2009):
If Anatomical Orientation Type (0010,2210) is absent or has a value of BIPED, the x-axis is increasing to the left hand side of the patient. The y-axis is increasing to the posterior side of the patient. The z-axis is increasing toward the head of the patient.
(we’ll ignore the quadupeds for now).
Three major cuts we have (or planes):
Transverse (AKA Axial) divides head from feet
Sagitall Cut - right between the eyes
Coronal Cut - the Filet
In a way it’s funny to call this the ‘patient-based’ coordinate system. ‘Doctor-based coordinate system’ is a better name.
DICOM defines a term: "Reference Coordinates System" or RCS. The RCS is a very intuitive coordinate system of the patient body: $X$ direction is from Right to Left. So if the patient is standing in front of you with the arm raised to the sides, then $X$ direction is from the right hand to the left hand.
$Y$ direction is from front to back or medical-wise from Anterior to Posterior so if the patient is standing in front of you so you see him/her from his/her left side, then $Y$ goes from your left to your right (confusing? look at the picture).
$Z$ direction goes from Feet to Head. At least this is simple to explain.
Now that we know the directions, there are letters assigned to the ends of each direction:
C.7.6.3.1.4 - Pixel Data
Pixel Data (7FE0,0010) for this image. The order of pixels sent for each image plane is left to right, top to bottom, i.e., the upper left pixel (labeled 1,1) is sent first followed by the remainder of row 1, followed by the first pixel of row 2 (labeled 2,1) then the remainder of row 2 and so on.
The resulting pixel array then has size (‘Rows’, ‘Columns’), with row-major storage (rows first, then columns). We’ll call this the DICOM pixel array. Thinking about it in the image plane (not with respect to the RCS), $X$ is vertical (rows) and increasing downward and $Y$ is horizontal (columns) increasing left to right.
Section 10.7.1.3: Pixel Spacing
The first value is the row spacing in mm, that is the spacing between the centers of adjacent rows, or vertical spacing. The second value is the column spacing in $mm$, that is the spacing between the centers of adjacent columns, or horizontal spacing.
See:
From section C.7.6.2.1.1 of the DICOM object definitions (2009):
The Image Position (0020,0032) specifies the $x$, $y$, and $z$ coordinates of the upper left hand corner of the image; it is the center of the first voxel transmitted. Image Orientation (0020,0037) specifies the direction cosines of the first row and the first column with respect to the patient. These Attributes shall be provide as a pair. Row value for the $x$, $y$, and $z$ axes respectively followed by the Column value for the $x$, $y$, and $z$ axes respectively.
From Section C.7.6.1.1.1 we see that the ‘positive row axis’ is left to right, and is the direction of the rows, given by the direction of last pixel in the first row from the first pixel in that row. Similarly the ‘positive column axis’ is top to bottom and is the direction of the columns, given by the direction of the last pixel in the first column from the first pixel in that column.
Let’s rephrase: the first three values of ‘Image Orientation Patient’ are the direction cosine for the ‘positive row axis’. That is, they express the direction change in ($x$, $y$, $z$), in the DICOM patient coordinate system (DPCS), as you move along the row. That is, as you move from one column to the next. That is, as the column array index changes. Similarly, the second triplet of values of ‘Image Orientation Patient’ ( img_ornt_pat[3:] in Python), are the direction cosine for the ‘positive column axis’, and express the direction you move, in the DPCS, as you move from row to row, and therefore as the row index changes.
Further down section C.7.6.2.1.1 (RCS below is the reference coordinate system - see DICOM object definitions section 3.17.1):
$$\begin{bmatrix} P_x \\ P_y \\ P_z \\ 1 \end{bmatrix} = \begin{bmatrix} X_x \Delta_i & Y_x \Delta_j & 0 & S_x \\ X_y \Delta_i & Y_y \Delta_j & 0 & S_y \\ X_z \Delta_i & Y_z \Delta_j & 0 & S_z \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} i \\ j \\ 0 \\ 1 \end{bmatrix} = M \begin{bmatrix} i \\ j \\ 0 \\ 1 \end{bmatrix} $$The Image Plane Attributes, in conjunction with the Pixel Spacing Attribute, describe the position and orientation of the image slices relative to the patient-based coordinate system. In each image frame the Image Position (Patient) (0020,0032) specifies the origin of the image with respect to the patient-based coordinate system. RCS and the Image Orientation (Patient) (0020,0037) attribute values specify the orientation of the image frame rows and columns. The mapping of pixel location ($i$, $j$) to the RCS is calculated as follows:
Where:
We stop to ask ourselves, what does DICOM mean by voxel $(i, j)$?
Isn’t that obvious? Oh dear, no it isn’t. See the DICOM voxel to patient coordinate system mapping formula above. In particular, you’ll see that thinking in terms of the traditional $X$-$Y$ graph:
That is, if we have the DICOM pixel data as defined above, and we call that pixel_array , then voxel ($i$, $j$) in the notation above is given by pixel_array[$j$, $i$]. Essentially we need to switch $i$ and $j$
What does this mean? It means that, if we want to apply the formula above to array indices in pixel_array, we first have to apply a column / row flip to the indices. Say $M_{pixelArr}$ is the affine to go from array indices in pixel_array to $mm$ in the DPCS. Then, given $M$ above: $$ M_{pixelArr} = M \begin{bmatrix} 0 & 1 & 0 & 0 \\ 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} $$
The (i, j), columns, rows in DICOM, columns, rows in DICOM) is rather confusing, so we’re going to rephrase the affine mapping; we’ll use $r$ for the row index (instead of $j$ above), and $c$ for the column index (instead of $i$).
Next we define a flipped version of ‘ImageOrientationPatient’, $F$, that has flipped columns. Thus if the vector of 6 values in ‘ImageOrientationPatient’ are $(i_1 ... i_6)$, then: $$ F = \begin{bmatrix} i_4 & i_1 \ i_5 & i_2 \ i_6 & i_3 \
\begin{bmatrix} F_{11} & F_{12} \\ F_{21} & F_{22} \\ F_{31} & F_{32} \\ \end{bmatrix} $$
Now the first column of $F$ contains what the DICOM docs call the ‘column ($Y$) direction cosine’, and second column contains the ‘row ($X)$ direction cosine’. We prefer to think of these as (respectively) the row index direction cosine and the column index direction cosine.
Now we can rephrase the DICOM affine mapping with:
$$\begin{bmatrix} P_x \\ P_y \\ P_z \\ 1 \end{bmatrix} = \begin{bmatrix} F_{11} \Delta r & F_{12} \Delta c & 0 & S_x \\ F_{21} \Delta r & F_{22} \Delta c & 0 & S_y \\ F_{31} \Delta r & F_{32} \Delta c & 0 & S_z \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} r \\ c \\ 0 \\ 1 \end{bmatrix} = A \begin{bmatrix} r \\ c \\ 0 \\ 1 \end{bmatrix} $$Where:
For later convenience we also define values useful for 3D volumes:
Let us say, we have a single DICOM file, or a list of DICOM files that we believe to be a set of slices from the same volume. We’ll call the first the single slice case, and the second, multi slice.
In the multi slice case, we can assume that the ‘ImageOrientationPatient’ field is the same for all the slices.
We want to get the affine transformation matrix A that maps from voxel coordinates in the DICOM file(s), to $mm$ in the DICOM patient coordinate system.
By voxel coordinates, we mean coordinates of form ($r$,$c$,$s$) - the row, column and slice indices - as for the DICOM affine formula.
In the single slice case, the voxel coordinates are just the indices into the pixel array, with the third (slice) coordinate always being 0.
In the multi-slice case, we have arranged the slices in ascending or descending order, where slice numbers range from $0$ to $N−1$ - where $N$ is the number of slices - and the slice coordinate is a number on this scale.
We know, from DICOM affine formula, that the first, second and fourth columns in $A$ are given directly by the (flipped) ‘ImageOrientationPatient’, ‘PixelSpacing’ and ‘ImagePositionPatient’ field of the first (or only) slice.
Our job then is to fill the first three rows of the third column of $A$. Let’s call this the vector $\textbf{k}$ with values $k_1$, $k_2$, $k_3$.
See also the definitions in DICOM affine formula. In addition
For the single slice case we just fill $\textbf{k}$ with $\textbf{n} \cdot \Delta s $ - on the basis that the $Z$ dimension should be right-handed orthogonal to the $X$ and $Y$ directions.
For the multi-slice case, we can fill in $k$ by using the information from $T^N$, because $T^N$ is the translation needed to take the first voxel in the last ( $slice index = N−1$ ) slice to $mm$ space. So:
$$ \left(T^N_1\right) = A \begin{bmatrix} 0 \\ 0 \\ -1 + N \\ 1 \end{bmatrix} $$From this it follows that: $$\left\{k_1 : \frac{T^1_1 - T^N_1}{1-N},\: k_2 : \frac{T^1_2 - T^N_2}{1-N} ,\: k_3 : \frac{T^1_3 - T^N_3}{1-N}\right\}$$
and therefore:
See derivations/spm_dicom_orient.py for the derivations and some explanations
We may have the problem (see e.g. Sorting files into volumes) of trying to sort a set of slices into anatomical order. For this we want to use the orientation information to tell us where the slices are in space, and therefore, what order they should have.
To do this sorting, we need something that is proportional, plus a constant, to the voxel coordinate for the slice (the value for the slice index).
Our DICOM might have the ‘SliceLocation’ field (0020,1041). ‘SliceLocation’ seems to be proportianal to slice location, at least for some GE and Philips DICOMs I was looking at. But, there is a more reliable way (that doesn’t depend on this field), and uses only the very standard ‘ImageOrientationPatient’ and ‘ImagePositionPatient’ fields.
Consider the case where we have a set of slices, of unknown order, from the same volume.
Now let us say we have one of these slices - slice $i$. We have the affine for this slice from the calculations above, for a single slice ($A_{single}$).
Now let’s say we have another slice $j$ from the same volume. It will have the same affine, except that the ‘ImagePositionPatient’ field will change to reflect the different position of this slice in space. Let us say that there a translation of $d$ slices between $i$ and $j$. If $A_i$ ($A$ for slice $i$) is $A_{single}$ then $A_j$ for $j$ is given by:
$$ A_j = A_{single} \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & d \\ 0 & 0 & 0 & 1 \\ \end{bmatrix} $$and ‘ImagePositionPatient’ for $j$ is: $$ T^j = \begin{bmatrix} T^1_1 + \Delta s d n_1 \\ T^1_2 + \Delta s d n_2 \\ T^1_3 + \Delta s d n_3 \end{bmatrix} $$
Remember that the third column of $A$ gives the vector resulting from a unit change in the slice voxel coordinate. So, the ‘ImagePositionPatient’ of slice - say slice $j$ - can be thought of the addition of two vectors $T^j = \textbf{a + b}$, where $\textbf{a}$ is the position of the first voxel in some slice (here slice 1, therefore $\textbf{a} = T^1$) and $\textbf{b}$ is $d$ times the third column of $A$. Obviously $d$ can be negative or positive. This leads to various ways of recovering something that is proportional to $d$ plus a constant. The algorithm suggested in this ITK post on ordering slices - and the one used by SPM - is to take the inner product of $T^j$ with the unit vector component of third column of $A_j$ - in the descriptions here, this is the vector $\textbf{n}$:
$$ T^j\cdot \textbf{c} = \left( T^1_1 n_1 + T^1_2 n_2 + T^1_3 n_3 + \Delta s d n^2_1 \Delta s d n^2_2 + \Delta s d n^2_3 \right) $$This is the distance of ‘ImagePositionPatient’ along the slice direction cosine.
The unknown $T^1$ terms pool into a constant, and the operation has the neat feature that, because the $n^2_{123}$ terms, by definition, sum to 1, the whole can be expressed as $\gamma + \Delta s d $ - i.e. it is equal to the slice voxel size ( $\Delta s $) multiplied by $d$, plus a constant.
Again, see derivations/spm_dicom_orient.py for the derivations.
In [1]:
""" Symbolic versions of the DICOM orientation mathemeatics.
Notes on the SPM orientation machinery.
There are symbolic versions of the code in ``spm_dicom_convert``,
``write_volume`` subfunction, around line 509 in the version I have
(SPM8, late 2009 vintage).
"""
import numpy as np
import sympy as sym
from sympy.abc import x
from IPython.display import display, Math, Latex
import IPython.display as disp
from __future__ import division
#from sympy.interactive import printing
"""after running this, all math output
will be formatted nicely for the duration
of the notebook session
"""
sym.init_printing() #use_latex='mathjax') #use_latex=True)
#sym.init_session(quiet=True)
#%load_ext sympyprinting
#%pylab inline
#import matplotlib.pyplot as plt
#%matplotlib inline
In [2]:
# The code below is general (independent of SPMs code)
def numbered_matrix(nrows, ncols, symbol_prefix):
return sym.Matrix(nrows, ncols, lambda i, j: sym.Symbol(
symbol_prefix + '_{%d%d}' % (i+1, j+1)))
def numbered_vector(nrows, symbol_prefix):
return sym.Matrix(nrows, 1, lambda i, j: sym.Symbol(
symbol_prefix + '_{%d}' % (i+1)))
# Dump out the formulae here to latex for the RST docs
def my_latex(expr):
S = sym.latex(expr)
return S[1:-1]
def spm_full_matrix(x2, y2, to_inv, inv_lhs):
rhs = to_inv[:,:]
rhs[:,1] = x2
lhs = inv_lhs[:,:]
lhs[:,1] = y2
return lhs * rhs.inv()
def disp_eqn(expr):
disp.display_latex(expr)
In [3]:
# premultiplication matrix to go from 0 based to 1 based indexing
one_based = sym.eye(4)
one_based[:3,3] = (1,1,1)
disp_eqn(one_based)
# premult for swapping row and column indices
row_col_swap = sym.eye(4)
row_col_swap[:,0] = sym.eye(4)[:,1]
row_col_swap[:,1] = sym.eye(4)[:,0]
disp_eqn(row_col_swap)
In [4]:
# various worming matrices
orient_pat = numbered_matrix(3, 2, 'F')
disp_eqn(orient_pat)
In [5]:
orient_cross = numbered_vector(3, 'n')
disp_eqn(orient_cross)
In [6]:
missing_r_col = numbered_vector(3, 'k')
disp_eqn(missing_r_col)
In [7]:
pos_pat_0 = numbered_vector(3, 'T^1')
disp_eqn(pos_pat_0)
In [8]:
pos_pat_N = numbered_vector(3, 'T^N')
pixel_spacing = sym.symbols(('\Delta{r}', '\Delta{c}'))
NZ = sym.Symbol('N')
slice_thickness = sym.Symbol('\Delta{s}')
disp_eqn(pos_pat_N)
disp_eqn(pixel_spacing)
disp_eqn(NZ)
disp_eqn(slice_thickness)
In [9]:
R3 = orient_pat * np.diag(pixel_spacing)
#disp_eqn(R3)
R = sym.zeros(4,2)
R[:3,:] = R3
#disp_eqn(R3)
disp_eqn(R)
In [10]:
# The following is specific to the SPM algorithm.
x1 = sym.ones(4,1)
y1 = sym.ones(4,1)
y1[:3,:] = pos_pat_0
disp_eqn(x1)
disp_eqn(y1)
In [11]:
to_inv = sym.zeros(4,4)
to_inv[:,0] = x1
to_inv[:,1] = sym.symbols('a, b, c, d')
to_inv[0,2] = 1
to_inv[1,3] = 1
disp_eqn(to_inv)
In [12]:
inv_lhs = sym.zeros(4,4)
inv_lhs[:,0] = y1
inv_lhs[:,1] = sym.symbols('e, f, g, h')
inv_lhs[:,2:] = R
disp_eqn(inv_lhs)
In [13]:
# single slice case
orient = sym.zeros(3,3)
orient[:3,:2] = orient_pat
orient[:,2] = orient_cross
disp_eqn(orient)
In [14]:
x2_ss = sym.Matrix((0,0,1,0))
y2_ss = sym.zeros(4,1)
y2_ss[:3,:] = orient * sym.Matrix((0,0,slice_thickness))
A_ss = spm_full_matrix(x2_ss, y2_ss, to_inv, inv_lhs)
disp_eqn(x2_ss)
disp_eqn(y2_ss)
disp_eqn(A_ss)
In [15]:
# many slice case
x2_ms = sym.Matrix((1,1,NZ,1))
y2_ms = sym.ones(4,1)
y2_ms[:3,:] = pos_pat_N
A_ms = spm_full_matrix(x2_ms, y2_ms, to_inv, inv_lhs)
disp_eqn(A_ms)
A_ms.simplify()
disp_eqn(A_ms)
# End of SPM algorithm
In [16]:
# Rather simpler derivation from DICOM affine formulae - see
# dicom_orientation.rst
# single slice case
single_aff = sym.eye(4)
rot = orient
disp_eqn(rot)
disp_eqn(pixel_spacing)
disp_eqn(slice_thickness)
rot_scale = rot * sym.diag(pixel_spacing[0], pixel_spacing[1],
[slice_thickness])
disp_eqn(rot_scale)
single_aff[:3,:3] = rot_scale
single_aff[:3,3] = pos_pat_0
disp_eqn(single_aff)
In [17]:
# For multi-slice case, we have the start and the end slice position
# patient. This gives us the third column of the affine, because,
# ``pat_pos_N = aff * [[0,0,ZN-1,1]].T
multi_aff = sym.eye(4)
multi_aff[:3,:2] = R3
trans_z_N = sym.Matrix((0,0, NZ-1, 1))
multi_aff[:3, 2] = missing_r_col
multi_aff[:3, 3] = pos_pat_0
est_pos_pat_N = multi_aff * trans_z_N
eqns = tuple(est_pos_pat_N[:3,0] - pos_pat_N)
solved = sym.solve(eqns, tuple(missing_r_col))
multi_aff_solved = multi_aff[:,:]
multi_aff_solved[:3,2] = multi_aff_solved[:3,2].subs(solved)
disp_eqn(multi_aff_solved)
In [21]:
#print(sym.latex(multi_aff_solved))
showing that the preivous equation is equivalent to the one shown earlier $$ \left[\begin{matrix}F_{{11}} \Delta{r} & F_{{12}} \Delta{c} & \frac{T^{1}_{{1}} - T^{N}_{{1}}}{1 - N} & T^{1}_{{1}}\\F_{{21}} \Delta{r} & F_{{22}} \Delta{c} & \frac{T^{1}_{{2}} - T^{N}_{{2}}}{1 - N} & T^{1}_{{2}}\\F_{{31}} \Delta{r} & F_{{32}} \Delta{c} & \frac{ T^{1}_{{3}} - T^{N}_{{3}}}{1 - N} & T^{1}_{{3}}\\0 & 0 & 0 & 1\end{matrix}\right]$$
In [23]:
# Check that SPM gave us the same result
A_ms_0based = A_ms * one_based
A_ms_0based.simplify()
A_ss_0based = A_ss * one_based
A_ss_0based.simplify()
disp_eqn(A_ms_0based)
disp_eqn(A_ss_0based)
assert single_aff == A_ss_0based
assert multi_aff_solved == A_ms_0based
In [25]:
# Now, trying to work out Z from slice affines
A_i = single_aff
nz_trans = sym.eye(4)
NZT = sym.Symbol('d')
nz_trans[2,3] = NZT
A_j = A_i * nz_trans
IPP_i = A_i[:3,3]
IPP_j = A_j[:3,3]
disp_eqn(IPP_i)
disp_eqn(IPP_j)
In [26]:
# SPM does it with the inner product of the vectors
spm_z = IPP_j.T * orient_cross
spm_z.simplify()
disp_eqn(spm_z)
In [29]:
# We can also do it with a sum and division, but then we'd get undefined
# behavior when orient_cross sums to zero.
ipp_sum_div = sum(IPP_j) / sum(orient_cross)
ipp_sum_div = sym.simplify(ipp_sum_div)
disp_eqn(ipp_sum_div)
In [37]:
str_latex = ('Latex stuff\n' +
' R = ' + my_latex(to_inv) + '\n\n' +
' L = ' + my_latex(inv_lhs) + '\n\n' +
' 0B = ' + my_latex(one_based) + '\n\n' +
' ' + my_latex(solved) + '\n\n' +
' A_{multi} = ' + my_latex(multi_aff_solved) + '\n\n' +
' A_{single} = ' + my_latex(single_aff) + '\n\n' +
r' \left(\begin{smallmatrix}T^N\\1\end{smallmatrix}\right) = A ' +
my_latex(trans_z_N) + '\n\n' +
' A_j = A_{single} ' + my_latex(nz_trans) + '\n\n' +
' T^j = ' + my_latex(IPP_j) + '\n\n' +
' T^j \cdot \mathbf{c} = ' + my_latex(spm_z) + '\n')
In [38]:
print(str_latex)
In [ ]: